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ABSTRACT 

We propose a rigorous, non-perturbative, Bayesian framework which enables one 
jointly to test Gaussianity and estimate the power spectrum of CMB anisotropics. 
It makes use of the Hilbert space of an harmonic oscillator to set up an exact likeli- 
hood function, dependent on the power spectrum and on a set of parameters ai , which 
are zero for Gaussian processes. The latter can be expressed as series of cumulants; 
indeed they perturbatively reduce to cumulants. However they have the advantage 
that their variation is essentially unconstrained. Any truncation (i.e.: finite set of ai) 
therefore still produces a proper distribution - something which cannot be said of 
the only other such tool on offer, the Edgeworth expansion. We apply our method to 
Very Small Array (VSA) simulations based on signal Gaussianity, showing that our 
algorithm is indeed not biased. 



1 INTRODUCTION 



With recent dramatic improvements in observations (Bernardis et al. 2000; Hanany et al. 2000) the cosmic microwave 
background (CMB) temperature anisotropics have become a gold mine of information in cosmology. In particular, the 
power spectrum Ci of CMB fluctuations has become a popular meeting ground for theorists and observers. On one 
hand measuring the provides a concrete target for improving experimental strategies. On another, well-founded 
theories predict that a system of peaks should be observed in the Ci , with a wealth of information encoded in their 



positions and heights. The recent unambiguous detection of the first of these peaks (Bernardis et al. 2000; Hanany 



et al. 2000: Miller et al. 1999) has caused a great deal of excitement (e.g. (Melchiorri et al. 2000 
2000| ; [Contaldi 2000D ) 



Dodelson & Knox 



Many reputable theories also predict that the temperature anisotropics form a Gaussian random field (even 
though there are notable exceptions (Salopek 199S; Linde & Mukhanov 1997; Peebles 1998; pontaldi. Bean & 



Magueijo 1999)). That being the case, the power spectrum does indeed contain all the relevant information char 
acterizing the random process under study. But how can we test the hypothesis of Gaussianity? Unfortunately the 
formalism used for Gaussianity tests is far less developed than its counterpart for power spectrum estimation. Even 
though Bayesianism has become the norm in current Ce estimation, tests of Gaussianity have all, but one ( pontaldi 
et al. 2000), been conducted in a frequentist vein. Most of these tests have revealed consistency with Gaussian- 



ity (e. g.(Kogut et al. 1996)), but claims for non-Gaussianity detection in COBE-DMR 4-y ear maps have also been 
made(Fbrreira, Magueijo fc Gorski 199§| ; Pando, Valls-Gabaud fc Fang 1998 ; Magueijo 2000). Of these, one (Ferreira. 



Magueijo fc Gorski 1998|) has been found to be due to an experimental syste matic (Banday, Zaroubi fc Gorski 1999 ) 
and another (Pando, Valls-Gabaud fc Fang 1998 ) to an error of method ( Barreiro et al. 2000 ). The third claim 



( Magueijo 200C ), however, remains in place. 

In this paper we set up a proper Bayesian framework with which to tackle the issue of non-Gaussianity. All 
Bayesian algorithms used so far have assumed that the signal, noise, and even galactic emissions are Gaussian 
random processes; and have targeted the Cg only. If the assumption of the signal's Gaussianity is to be dropped then 
one must propose an alternative form for the likelihood. Here the difficulties begin, since a workable non-Gaussian 
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likelihood cannot be easily defined. Non-Gaussianity spans an infinite dimensional space, and truncations into a 
finite number of degrees of freedom are usually inconsistent. 

The only attempt made so far is published in (Contaldi et al. 2000). In that work the Edgeworth expansion 
(Kendall & Stuart 1977) was used to produce a Bayesian joint estimate of the and the skewness of the COBE-DMR 
4-year maps. This approach is at best very approximate. Cumulants of distributions are subject to very complicated 
constraints (Kendall & Stuart 1977), and setting them all to zero apart from one (in (Contaldi et al. 200C) the 
skewness) is inconsistent with the basic properties of a distribution, such as positive definiteness or normalization. 
Furthermore the Edgeworth expansion is an asymptotic series, which is never a distribution when truncated. 

In this paper we bypass these difficulties by deriving the general form of the likelihood from the Hilbert space of a 
linear harmonic oscillator. We recall (Merzbacher 1970) that the ground state has a Gaussian wave- function, while all 
excited states have wave-functions which multiply a Gaussian by a Hermite polynomial. The space of all distributions 
can then be spanned by the amplitudes, a„, of the various energy eigenstates, with a general distribution taking the 
form of a Gaussian times the square of a (possibly finite) series of Hermite polynomials. 

Such a generic distribution has a remarkable similarity with the Edgeworth expansion, which takes the form of 
a Gaussian multiplying an infinite series of Hermite polynomials with coefficients which are themselves polynomials 
in the cumulants of the distribution. Closer comparison of the two expressions reveals that indeed the amplitudes 
a„ can be written as series of cumulants( |Contaldi, Bean fc Magueijo 199"e ); these are the combinations of cumulants 
which can be varied independently. In particular these are the combinations which can be independently set to 
zero without mathematical inconsistency. Furthermore, perturbatively (that is when the cumulants are "small" in 
a suitable sense), the amplitudes a„ are proportional to the n*^ order cumulant. In some sense the a„ generalize 
cumulants to non-perturbative situations. 

We thus arrive at a well-defined mathematical framework for conducting Bayesian tests of Gaussianity, which 
jointly produces Ci estimates. Its interest is twofold. Firstly there is the obvious interest in finding out whether the 
CMB fluctuations are Gaussian or not. Secondly there is the issue of whether Ce estimates themselves may vary if 
non-Gaussian degrees of freedom are allowed into the likelihood. In this paper we describe this formalism (Section 
and apply it to VSA simulations (Section pending actual data. 

It is interesting to note that the are more than a mathematical device, and have a physical interest of their 



own within the framework of the standard inflationary scenario(Guth 1981; Lindc 1982; Albrccht & Steinhardt 1982 



Linde 1983). Standard inflationary fluctuations are Gaussian because the infiaton's fiuctuations satisfy harmonic 
oscillator equations, and are assumed to be in the ground state. The latter is an assumption which relics loosely 



on the boundary conditions imposed in quantum cosmology (Hartle & Hawking 1983; Halliwell & Hawking 1985 



Vilcnkin 1982), and needs not be correct(Contaldi, Bean & Magueijo 1999). A non-trivial wave-function for the 
infiaton's fiuctuations manifests itself in non-Gaussian density fiuctuations, even if we do not depart from single-field, 
slow-roll, infiation. Quantifying their non-Gaussianity by means of the a„ offers a direct bridge to the wave-function 
of the infiaton's fiuctuations. 

Hence, if we take it for granted that infiation is realized in its simplest form, and is triggered by quantum 
cosmology, then measuring the parameters a„ amounts to mapping the wave-function of the Universe as it emerged 
out of a quantum epoch. 



2 SIGNAL-TO-NOISE EIGENMODES 



The non-Gaussian likelihood formalism we are about to present works most simply when applied to a series of 
independent variables. We shall therefore combine our method with the technique of signal-to-noise eigenmodes, 
which we start by reviewing. 

The signal-to-noise technique is a special case of the Karhunen-Loeve method where the parameter dependence 
is linear (affine) (Bond 1995; Bond, Jaffe & Knox 1998; Tegmark, Taylor & Heavens 1997). Let us consider a general 
set of random variables 



(1) 



where Si is the signal component and Ui is the noise counterpart. Let us also assume that the signal and noise 
contributions are independent and each has zero mean. The quantity di could represent, for example, the observed 
temperature fiuctuation in the ith pixel of a CMB map, or alternatively the real or imaginary part of the amplitude 
of the ith coefficient in the Fourier expansion of the map. The covariance matrix {didj) is given simply by Cij = 



{Si 



— Si- 



Nij , where S and N are the signal and noise covariance matrices respectively. 



In the standard likelihood approach one estimates the parameters oi, 02, . . . , a„ of the probability distribution 
from which the Si are drawn by calculating the likelihood £(d|a) as a function of these parameters. The parameters 
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a usually enter the calculation through the signal covariance matrix S'(a). Let us suppose, in some part of the 
calculation, we are interested only in the parameter ai (say), with the remaining parameters fixed at particular 
values. In some (very common) cases, it proves possible to write the signal covariance matrix in the special form 

5(a) =5i(ai) + ^2(02,..., a„), (2) 

where Si depends only on the parameter of interest and S2 depends only on the (fixed) values of the other parameters 
and is thus a fixed matrix. Furthermore, it is often also the case that 5*1 is linear in ai and so can be written in the 
simple form 51(01) — aiU where the fixed matrix U — 6*1(1). Thus the total covariance matrix C can be written as 

C = S + N = aiU + S2 + N = aiU + V, (3) 

where U is the 'unit signal' covariance matrix and V is the 'generalised noise' matrix. In particular, we note that if 
the parameters a are the power spectrum coefficients (or averages of the C^'s in given ^-bins), it is always possible 
to write the covariance matrix in the form (^. 

Since U and V are both real symmetric matrices, they can be diagonalised simultaneously by a single similarity 
transformation. This is most easily achieved by solving the generalised eigenproblem [/x = Al^x. Let us denote the 
corresponding eigenvalues by A; and eigenvectors by e.^, which are normalised such that e\Vei = 1. If we now consider 
the new set of variables £,i — Gi ■ d, then it is straightforward to show that these new variables are uncorrelated for 
any value of the parameter ai, with a covariance matrix given by = (1 + ai\i)5ij. The are the so-called 

eigenmodes of signal-to noise (S/N); the modes with large eigenvalue are expected to be well-measured while modes 
with small eigenvalues are poorly-measured and do not contribute significantly to the likelihood. 

If the original data di were distributed as a multivariate Gaussian, the particular advantage of the S/N eigenmode 
basis is that the likelihood function for the parameter ai (with 02 , . . . , a„ held fixed) becomes a simple product of 
one-dimensional Gaussians, and can be written simply in terms of the eigenvalues as 

mai)^l[ 1_ e-^i^. (4) 

-^.-^ VSttvI + aiAi 

It is clear that this procedure can be repeated for each of the other parameters 02, . . . , a„, provided in each case 
the covariance matrix can be written in the form (|^) . The likelihood function can thus be evaluated very simply along 
the 'coordinate directions' in parameter space. Moreover, in the special case where the parameters oi, a2, . . . , a„ are 
mutually independent (or quasi-independent to a good approximation), the likelihood function factorises as 

/:(d|a) = C{d\ai)Cid\a2) ■ ■ ■ £(d|a„). (5) 

Thus, in this case, the above signal-to-noise eigenmode procedure can be repeated for each parameter in turn to 
evaluate the full likelihood function. 



3 AN EXACT, NON-PERTURBATIVE, NON-GAUSSIAN LIKELIHOOD 

Let X represent a general random variable, within a set of variables which are assumed to be independent. Let us build 
its distribution from the space of wave- functions which are energy eigenmodes of a simple harmonic oscillator. The 
following results may be found in any quantum mechanics book (e.g.( Merzbacher 1970| )) adopting the Schrodinger 
(rather than the Heisenberg) picture. We have that the general wave-function for x is given by: 

■^(a;) = ^ (x) (6) 

n 

where n labels the energy spectrum E = hu!{n + 1/2). The basis functions ipn take the form 



ijn{x)^CnH,A -=— e (7) 
V v2cro/ 

with normalisation fixing C„ as, 

Cn = X= r • (8) 

(2"n!V2^ao)^/2 

The only constraint upon the amplitudes a„ is: 

n 

This is a simple algebraic expression which can be eliminated explicitly by writing: 
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ao 



\ 



El 



(10) 



The quantity ctq is the variance associated with the (Gaussian) probability distribution for the ground state |-0o| 
We shaU work with Hermite polynomials Hn{x) defined as 



Hr^ix) = (-l)"e- — e-- 
and normalised as 

The most general probability density for the fluctuations in x is thus: 

2 



Cfi Hrt 



V2ao 



(11) 



(12) 



(13) 



The ground state (or "zero-point" ) fluctuations are Gaussian, but any admixture with other states will be 
reflected in a non-Gaussian distribution function. Accordingly we may use the amplitudes of these admixtures, a„, 
as non-Gaussianity indicators. Their obvious advantage is the rather trivial constraint (^), which can be ignored using 
(|lO|). It permits concentrating on a finite set of non-Gaussian degrees of freedom, without mathematical inconsistency. 

However there is another reason why the a„ are of mathematical interest: they reduce to cumulants k„ under 
certain assumptions. If we assume mild non-Gaussianity (which we define through the condition jaoP S> Ictip, for 
i > 1) then to first order in a^. 



Pix) 



2tt(To 



1 + 



n> 



^^(2«n!)V2 "VV2ao 



(14) 



where we have taken ao to have zero phase (so that to first order ao = 1)- Comparing (|lj) with the Edgeworth 
expansion (Kendall & Stuart 1977) we find a one-to-one correspondence between the amplitudes of the various 
energy eigenstates, and the combinations of cumulants appearing as coefficients in the Edgeworth expansion. The 
latter simplify enormously if we only keep first order terms, that is if we assume that quadratic and higher order 
terms in the cumulants are negligible. Then we find that 

K„ cx 3fJ(a„) (15) 

with a rather complicated proportionality constant (which is easy to work out case by case). Hence, to first order, 
the coherent contamination of the ground state by the n*'' energy eigenstate is signalled by a non- vanishing cumulant 
Kn- For instance the presence of the third energy eigenstate results in K3 oc 5R(a3) 7^ 0, and, to first order, zero higher 
order cumulants. 

The advantage of using the a„ is that they still work (i.e. they still lead to proper distributions) when the 
distribution is highly non-Gaussian. Any maximum likelihood estimate will necessarily have to probe regions of a^ 
beyond the perturbative regime, even for Gaussian realizations. In these regions setting all but a finite number of 
K„ to zero is inconsistent; but not for a finite number of a„. In the non-perturbative regime the a„ become rather 
complicated series of k„. However these series of cumulants may be varied, or set to zero, independently, and still 
lead to a distribution. Hence we should regard the a„ as non-perturbative generalizations of cumulants. 

More concretely, in the non-perturbative regime, we have 



P{x) = IV^I^ 



27r(To 



E 



H, 



( ^ 



(2»+J'z!j!)i/2 



H 



\\/2< 



Co 



We may recover the Edgeworth expansion by noting that. 



e-" H,{x)H,{x) = e- 
with, 

" [s — n)\{s ^ i)\{s — i)\ 



(16) 



(17) 



(18) 
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with 2s = n + i + j. One may derive (|l^) usi ng (p^) and the standard r esult for the integral over a product of 
three Hermite polynomials (formula 7.375.2 of ( Gradshteyn fc Ryzhik 199^ )). Thus we obtain the more complicated 
expression 



Comparing with the Edgeworth expansion leads to the rather complicated non-perturbative expression relating the 
amplitudes Ui and series of cumulants. 



4 APPLICATION TO VSA SIMULATIONS 



We have applied our method to simulated observations by the Very Small Array (VSA) interferometer. The VSA 
has been built by Cambridge and Jodrell Bank in the UK, and is sited at the Teide Observatory in Tenerifc. It 
has just become operational. The VSA is expected to give detailed maps of the CMB anisotropy with a sensitivity 
~ b^K and covering a range of angular scales from 10' to 2° for a frequency range of 28 and 38 GHz. It will have 14 
antennas and a 2-GHz bandwidth analogue correlator and uses the same technology as Cosmic Anisotropy Telescope 
(CAT). The VSA is able to observe in 'compact' and 'extended' modes, which are sensitive to different ^-ranges. In 
the compact mode, it is expected that VSA will recover the angular power spectrum in 10 spectral bins where each 
bin is centred respectively at ^ « 114, 211, 308, 404, 501, 598, 695, 792, 889, 986. The width of each bin M « 97 
corresponds to the 1/e diameter of the aperture function of the interferometer and represents the length scale on 
which the underlying Fourier modes of the sky are correlated by the instrument. Thus, the spectral bins have been 
chosen so that the power spectrum estimates in each bin are quasi-independent to a good approximation. 

An acco unt of a maximum-likelihood met hod for analysing interferometer observations of the CMB anisotropics 
is given in ( Hobson, Lasenby fc Jones 1995 ). In the standard likelihood analysis the parameters of interest are 
ctfc = {^"^ C I / {2it)) i^^Yi bin- Since these parameters are quasi-independent, the full likelihood function factorises as in 
(^ to a good approximation. The signal-to-noise eigenmode procedure outlined in section ^ is then applied to each 
factor, which, in the standard approach, reduces to the simple product of one-dimensional Gaussians given in (||). 

In our new approach, however, we subject the standard Gaussian likelihood algorithm to the following modifi- 
cations. Instead of assuming the simple Gaussian form for the probability distribution of each S/N eigenmode ^i, we 
instead consider the more general situation in which all a„ are set to zero, except for the real part of a^. The reason 
for this is that such a quantity reduces to the skewness in the perturbative regime. The imaginary part of aj, is only 
meaningful in the non-perturbative regime (and can be set to zero independently without inconsistency). Hence we 
are considering a likelihood of the form: 



P{x) 
with 



27r(To 



1 



CkO H 1=^1^, 



^2, 



(Jo 



n 2 



(20) 



(21) 



explicitly replaced in (|20|). 

The generalization of this distribution to the multidimensional case is trivial in the signal-to-noise eigenmode 
basis, since we can simply take the product of the individual one-dimensional distributions. Thus, when considering 
the power spectrum in the fcth spectral bin, we adopt the likelihood function 



-^(^10?=, as) = W 



2(l + <.fcAi) 



27r\/l + flfeA,; 



as 



(22) 



where ak is the average of value l'^Ci/{2Ti) in the fcth spectral bin. The could in principle depend on t, but for 
simplicity we have dropped this dependence. Notice that because the noise is Gaussian, and because of the principle 
of superposition in quantum mechanics, the a„ appearing in this formula are the ones pertaining to the signal alone. 

We have applied this method to a 30 x 12 hour simulated VSA observation of a Gaussian CMB realisation 
drawn from a standard inflationary model with rJcdm = 0.95, Ob = 0.05, JIa = 0, /i = 0.5, = 1 and no tensor 
contribution. In Fig. 1 we show the contour plots of the likelihood functions for the amplitude of the power spectrum 
€^Cf/(27r) and for each of the 10 spcctrals bins observed by VSA. The alignment of the contour axes with the 
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coordinate axes implies the reassuring result that there seems to be little correlation in each bin between the power 
spectrum estimate and the estimate of the non-Gaussian parameter. 

In Fig. 2 we plot the likelihood function for in each spectral bin after marginalization over the amplitude 
of the angular power spectrum. The results obtained indicate that the value of as scatter around as = 0, within 
the range implied by the width of the likelihood. The percentage of the population inside the contour intersecting 
the origin represents the confidence level for rejecting Gaussianity. All of these are within 1-2 sigma, indicating that 
this method is not biased. Notice that for those bins in which we failed to obtain a CMB detection (see Fig. 3) there 
seems to be a bias towards a peak at — 0, without the scatter expected from the width of the likelihood. 

In Fig. 3 we plot the likelihood functions for the amplitude of the power spectrum in each bin after marginalization 
over the parameter (solid line) . Superimposed are the corresponding conditional distributions for 03 = (dashed 
line). For this Gaussian CMB realisation, the distributions obtained are not significantly affected by the inclusion 
of the extra parameter 03 . The most noticeable effect is a slight variation of the position of the peak (particularly 
for bins 1 and 2) which is in agreement with Fig 1. Since in each bin the estimate of the power spectrum and the 
non-Gaussian parameter are weakly correlated, we see that the widths of the likelihood functions for the angular 
power spectrum are not significantly increased by including the ^3 parameter. 

In Fig. 4 we plot the joint likelihood for 0:3 obtained by multiplying the individual likelihoods for the 10 spectral 
bins, thus obtaining the overall estimate and a better constraint on a^. 

Finally, in Fig 5 we plot the distribution of the peak of the likelihood for obtained from Monte Carlo simu- 
lations. In each VSA simulation the CMB is a realization drawn from an inflationary model. The CMB fluctuations 
are thus Gaussian. The distribution peaks around a value of q;3 = confirming that our algorithm is indeed not 
biased. 



5 CONCLUSIONS 

In this paper we laid down the foundations for a rigorous Bayesian framework for testing non-Gaussianity, and jointly 
estimating the power spectrum (Section |^) . Our main achievement was to convert testing Gaussianity into a problem 
of Bayesian estimation. We defined a series of parameters a„, to be added to the power spectrum, such that the 
origin of the new space represents Gaussianity. These parameters are generalizations of cumulants. If all cumulants 
are very small, in a suitable sense, each of the new parameters is approximately proportional to a cumulant. If not, 
then the new parameters become series of powers of cumulants. They are desirable non-perturbative generalizations 
of cumulants because they are independent, ie: subject to essentially no constraints, unlike standard cumulants. 

With any dataset, one must then determine the contour of the likelihood intersecting the origin of the a„ space, 
after marginalization over the power spectrum. The percentage of the population inside this contour is the confidence 
level for rejecting Gaussianity. We found that for simulated VSA observations of a Gaussian CMB realisation this 
confidence level is always within 1-2 sigma. To assess if our algorith is unbiased one must produce simulated VSA 
observations of several Gaussian CMB realizations. We found that the distribution of the peak of the likelihood in 
as for a number of these realizations peaks around a value of 0:3 = showing that our method is indeed not biased. 
(Section |). 

The method we have proposed is completely general, and may be applied to any type of experiment, interfer- 
ometric or single-dish. In particular its application to COBE-DMR maps, closely mimicking the steps of ( |Contald"i 



200c), is straightforward. The only issue which may complicate the method is galactic foreground removal. In some 
experiments foregrounds away from the galactic plane may be ignored, by suitably choosing the frequency channels. 
In some cases, contaminations may be removed by subtracting off the correlated component, making use of templates 



(Kogut et al. 1996). In these cases there is no extra complication to our method. 



However in some cases (Hobson, Lasenby & Jones 1995) foreground subtraction is part of the maximum likelihood 



algorithm leading to CMB Ci estimates. In some of these cases it is assumed that Galactic foregrounds form a 
Gaussian random field. With our method we may now allow for non-Gaussian degrees of freedom to be applied to 
these emissions. Hence we should be able to improve significantly on these methods of foreground removal, as well 
as exploring signal non-Gaussianity. The detailed implementation of this project, as well as its application to VSA 
data, is the subject of a future publication. 

The formalism we have developed is also of assistance for generating realizations belonging to the most general 
ensemble parameterized by the an- In work in preparation we show how this can be done, and how the maximum 
likelihood method proposed in this paper may then differentiate between distinct distributions on the basis of single 
realizations. 
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Figure 1. Confidence contours for a simulated VSA observation of a Gaussian CMB realisation drawn from a standard inflationary 
model. In each spectral bin the fluctuations are parametrized by the average value of C; / (27r) and the generalised cumulant as . Contours 
are at 10,20,30, ...90 and 95 per cent confidence levels. 
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Figure 2. Marginalised likelihood functions for as in each spectral bin obtained from the simulated VSA observation. 
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Figure 3. Solid line: marginalized distribution for CMB fluctuations parameterized by PCi /{27r), for a simulated VSA observation; 
dashed line: conditional distribution for as = 0. 
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Figure 4. The joint likelihood for as obtained by multiplying the individual likelihoods in 0:3 for the 10 spectral bins. 
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Figure 5. The distribution of the peak of the likelihood in as for simulated VSA observations of several Gaussian CMB realizations. 



